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MAY 1, 1976 - JULY 31, 1970 
SONIC BOOH RESEARCH 


During this period we completed our numerical program for the sonic boom 
analysis including the asymmetric effect due to lift near the vertical 
plane of _. r i,imetry. Several numerical examples were computed. The results 
were presented at the third Aeroacoustic Conference in July (Two copies of 
AIAA Preprints No. 76-937 are attached). ' 

The program and a numerical example were delivered to the contractor. A 
description for the usage of the program is presented In NYU Report AA/76-1 1 
entitled "Numerical Program for Sonic Doom Analysis-Nonlinear With Asymmetric 
Correction Due to Lift, by Fanny Rung, (Two copies are attached). 


NUMERICAL PROGRAM TOR SONIC DOOM ANALYSIS* 

- NONLINEAR WITH ASYiiMl.TRlC CORRECTION DUE TO LIFT 


F, RUNG** 

DEPARTMENT OF APPLIED SCIENCE 
NtU YORK UNIVERSITY 
NEW YORK* NEW YORK 100Q3 
NYU/AA76-11 

* * * ‘k t. -k i: * k * •!: * *V * -k k •!: 

ABSTRACT 

A coiiiputer program for CDC 6600 Is developed for the nonlinear sonic boom 
analysis including the asymmetric effect of lift near the vertical plane 
of symmetry. The program is written in FORTRAN IV language. This program 
carries out the numerical integration of the nonlinear governing equations 
from the input data at a finite distance from the airplane configuration at 
a flight altitude to yield the pressure signitude at ground. The required 
input data and the format for the output are described. A complete program 
listing and a sample calculation are given in the Appendix.^ 


* This research is supported by NASA Grant 
*'■' Research Scientist, Department of Applied 

(1) Appendix will be forwarded upon request. 
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LIST OP SYMBOLS 


loco! speed of sound In ft/sec. 
speed of sound at: flight altitude 

outgoing and incoming characteristics In the plane » o. 
shock shape, x « r(r, i|i) 

32.2 .. ft/sec 2 
altitude in ft. 

reference length (&cy the body length) in ft. 
local Mach number 
flic, Mach number 

pressure increment from ambient value in lbs. /ft? 
radical uistance/L 
gas constant 
entropy 

S /R ~ entropy induced by shock wave/R^at t[< - o 
(S sh /R) # at * = ° 

axial velocity radial and circumferential components in ft,/ 

flight velocity - M a In ft. /sec. 

o o 

u # /U at * - o 

v^/U at tji •=• o 

w./U at , - o 

temperature in °Ranlcine 

axial coordinate (hor izontal)/L 

vertical caordinate/L 

f 1 ight altitude in ft. 


List oi Symbols (Coni'* d) 


a 

3 

Y 
0 

V 
v 


the; Iwlf cone angle of tin; equivalent body in degree 

shock angle with the flow direction ahead of the shock 

ratio of specific heat ” \.h 

flow inclination with x-oxls 

Mach angle 

Prandtl Meyer angle 

meridian angle, i|) e o is the vertical plane of symmetry 
below the airplane 


Subscripts 

» the undisturbed value along the same streamline 

x,r,i}> partial derivatives 



outlim: or mr humfrical mmv, 


Tho basic equations and the justifu <ns for the approximations in the 

asymmetric effects for this numerical program were described in ref. 1, ?, 

** ** *. •* 

The numerical analysis deals with nine variables, M, o, r w , S^ , u, v, w, y fj 

r* 

and S s | as functions of x, r in the vertical plane of symmetry ^ - 0. < ,e 

first four variables are the variables in the quasi -symmetrical analysis 
(ref, 1) and the remaining five variables, which are the circumferential derive 
tives, u /U , v, ,/U U w,, (y ),./! and (S .) /\\ at f *• o, represent i ng tint 
asymmetric variables (ref. 2). The nine governing equations for these two 

A 

sets of variables are coupled weakly through one variable w. 

The iterative procedure for the numerical integration of those nine equations 

A 

begins with the assignment of the value w. The four equations for the four 

r* 

variables M> o, r ro , are determined in the same manner as in the axi -sym- 
metric calculations (ref, 1). The location of the new grid point x> r and 
the values M, Q are given by the characteristic equations with a curved C 
line segment* The values for r^ and are obtained by integrations along 

the stream line. We will then proceed to compute the asymmetric variables. 

+* . ** •[- „» 
u and v will be determined by the integration of the equations along C and c 

A M M* 

w» y ra and $ , are obtained by the integration of the equations along the 
streapi line. With this new w we continue to the next cycle of iteration until 
convergence. 

Similar iteration scheme is employed for the shock point as described in 
ref. 2. 

Since the asymmetric computations can be considered as the determination of 

A 

the value of w in the quasi -asymmetric computations. The general scheme of 

the computation from the initial data to the final pressure signature u ; rouni 

- 1 - 


will bo the same as that for the axi -symmetric computations; 


The calculations proceed along successions of C" linos with index I. Along 
each c" line, it will compute the new bow shock point and then the chnrae- 

j. 

teristic points with index J until the last C line with index ]V or the 

IllU A 

C line from a input data point. Then we proceed to the next c line with I 
increased by 1. 

When two C 1 * lines cross over, an embedded shock is formed. ' The integer KS is 
assigned in the order of shocks formed. For each embedded shock, its location 
J along the C~ line is denoted as J ( KS ) . In each sweep along a C” line, when J 
equals J ( KS ) we call the subroutine for embedded shock point (EMSHOCK) instead 
of C + - C“ (CPCM) subroutines. 

Whenever a shock front or a C line cuts across the ground (r “ Y/i.) the values 
on the ground level are obtained by linear interpolation. The pressure signa- 
ture includes the reflection coefficient 1.8. When the grid point from the 
last C line is below ground, the computation ends. 

The program contains four main subroutines; ( i } c 1 and bow shock [CPSH], 

(ii) C 1 and C~ [CPCM], (iii) formation of shock [FSHOCK] and (iv) embedded 

shock [EMSHOCK] . It does not have a subroutine for the intersection of two 

shocks therefore the program will stop when two shocks do intersect. 

1 

The input definition will be discussed in Table I, the format for input data 
appears in Table II, and the description of output will be in Table III. 

We will discuss briefly how do we control the step sizes and prepare the in- 
put data. 



COMTRUL OP STIP SIZES 


The step size is controlled by Inequality • > in the radio! distance vjr.cn 

two adjacent bow shock points? r*r<Cj for r' £ ?i 

r - r < Cj r" for < r" < r ? 

_** „ « 
and \ - r < C ? , for r ?. < r 

In this program we have f* =1 c : “0,04 r 2 « 100 Cj, - 4. Therefore 

the upper bound for the stop sizes increases linearly from 0,04 at r - fi" 1 

to 4 at r - r 2 = 100 and then rental ns equal to 4 all the way to the ground. 

A change in step size and in the rate of increment can be made by changes 
in the values of Cj r 2 and C;> provide that it is continuous* ,i.c., 

Ci r ? * C 2 . 

When we want to use different control functions wo should change the control 

*{* 

equations in the main program and in the subroutine for C and shock (CKSil), 


PREPARATION or INPUT DATA 


The input data along r * r Q will be obtained from experimental data or 
from the full three dimensional analysis near the airplane configuration. 

M « •< A fw 

At each point we input x> H, o, r u , S h , u, v* w, y* and S^, for the 
bow shock we input in addition the shock angle f>. It should be observed 
that not all the input data are independent of each other. For example: 
at bow shock* r^ - r while S j, P> and M ore. related to o and f and S ^ 
is re laved to y m > When tiie input data are read from the data cards, we get 

IUCAL « 0 

If we intend to compute or modify some of the input data by some equations 
we set 

IUCAL « 1 

and make the modifications from Statement no. 107 to Statement No. 3505. 

In our sample calculation, the input data will first be prepared from an 
axi -symmetric computation of an equivalent symmetric body to r - ? , There- 
fore we have u = v = w = S gh - 0, Since y = Y/L - r cos i> we have y^ “ r tl . 

If we sot IUCAL = 0, the pressure signature will be that of a symmetric body 

with asymmetric effect due to the two dimensional atmosphere layer only and 
the difference from the pure ax i -symmetric calculation is very small. 

We set IUCAL - 1, as in the sample calculation and then proceed to compute 

•* ** A +* *» 

the asymmetric terms u, v, w, S ^ and y ro from the linearized theory of an 

assigned lilt -i 'ribution (Ref. 2). Since the leading edge of wing is 
so located, that the characteristic line hits r = r at x = 6.85 in the plane 
ip » o lying behind the bow shock, the bow shock is still symmetric with S r ^ ~ 0 

rw f.. A ^ 

and y^ - r ot , Only the input data of u, v and w for x > 6.85 are changed 



by the formulas for the special lift variation in Fig. B oi ref. 2. for <» 
different approximation theory or a different lift variation tlio'.e equations 
should bo revised, (Statement flu. 170-3W5), 


Finally, we want to point out that we have included the data beyond the body 
in order to yield the tail shock at ground. In the sample calculation we 
use five body lengths. An estimate of the length on the safe side can be 
obtained from the length required in the real atmospheric program based on 
Whitman's theory. 
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1. Terri, A. j Siclmi, M. «tnd Ting* L : Sonic Deem Analysts for High 
Altitude flight at High Mach Number, AIAA Paper 73-1034, Progress 
in Astronautics and Aeronautics Vol. 38 pp. 301-320, AIAA, [tv; Y ork , 
1970 

2. Terri, A.; Ting, L. and to, R. Nonlinear Sonic Doom Analysis Including 
the Asymmetric Effects. AIAA Paper No. 76 - Aero acoustlc Conference 
Stanford, Calif,, July 26-28, 1976. 
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XXI 10 -™-— - 


TEST R 


TABLE 1 


INPUT DEnumORS 


flight altitude in ft, 

reference length (say the body length) in ft. 
flight Mach number 

radial distance of initial data lino/L 


TEST IM } not used In this program 


XXXXX1 


" C1 1 


XXXXX2 , > constants in the control of step size 


— 


« 0 


IUCAL 


isrop 


NUMBE 


IUMBE 


all input read from data cards 

some input data compute from formulas for each 
special cases 

maximum number of c line, for I stop use 
a large number* say 1000 to go to the ground. 

for detail information along I*th c* v line, sot ISTQP ~1* 

y , altitude in ft. 

Is, temperature in °R 

number of HH table (input card 5). In this program 
I use number = 63. In case the altitude is greater 
than 120,000 ft. We can increase this table by changing 
the value of NUMBE. 

Prandtl Mayer angle 

M> Mach number 

Number of entries in VV table (input card Ci). In this 
program I use IUH.'jf = 172, In case the Mad* number is 
greater than 3.5, the table has to be increased and. so 
does the value of IUMBE. 


number of embedded shock 
3 » of the (‘itibodded shock 
3 * of Dow shock 
axial coordinate (horizontal) /L 
vortical coordinate /L. 
flow inclination with x-axis. 
local Mach number 

M 

value of r for a stream lino fnr upstream 
entropy due to shock wave (s)/K 
u . . /U cat <f' = 0 

ijujr tn 

c 0 

w „, / u . It ♦ ■ 0 
y ^./ 1 - ot *" 0 

< s ,hV* st «*“ 0 


used in the special example for creation of 
asymmetric data with IUCAL = 1. 





table n * 


C ard 1 

Card 2 

Card 3 
Card 4 
Card S 

Card 6 

Card 7 

Card 8 

Card 9 
Ca rd 10 


INPUT DATA FORMAT 
(4F 10,5) 

H 0 » XL, XXMQ s CRD 

(3F 10,5) (not ^eing USG£ i J put a 

TESTR, TESTIM, TESTX blank card) 


(/IF 10,5) 



XXXXXX1, XXX XX?., XR1» XR2 



(415) 




IUCAL, ISTOP 



(2F 10,5) 

HH(1) 

1 

1 

TT(l) 

1 

1 



1 

HH(NUMBE) 

1 

TT(NUMBE) 



(10X, 2E 

20.8) 



W(l) 

ZMMf'U 1 ) 



W(2) 

1 

1 

ZMMM(2) 

I 

1 



t 

VV(IUMBE) 

t 

ZMMH( IUMBE) 



(5X, 15, E20. 10) 

MAX BOV , DELTA (not being used, 

put a blank card) 


(15, E15.8) 
KS, 

BETT 



(4E2Q.9) 

BIGXTEM, 

' EPSTEM,. SBJ.RKK, 

JMs’lR 


(3E 20.9/ 10X, 3E20.9) 


XX(1) 

THETHE(l) 

XMXM(l), RlflRIN(l) , 

SRSR(l) , BET 

UU ( 1 ) 

, VE(1) 

WV/ (1) , YY ( 1 ) , 

SS (D 

XX(2) 

THETHE(2), 

XMXM (2 ) , RINRIN (2) » 

SR5R (2 ) 

UU (2) 

VE(2) , 

WW(2) Y V ( 2 ) , 

SS (2) 


TABLE IN 



s 


I 



I 

[ 

( 


I 

l 

i 


! 

! 

1 

I' 



The output for this program is divided into two parts. The first port con- 
tains a tabulation of the input parameters to the program together with 
throo input tables, JIH vs. TT; VV v.s. ZMMM; and the input table of x,0 M , 
r co $ | » wltli the given radius (CRB). The second part first contains with 
table of J, u , v, w» y w , which will be the values of the input data 
if IUCAL * 0. in case IUCAL y* 0, the values will be recalculated (sec 
remarks In the preparation of input data) and the recalculated values of u, 
v, w, y M , S s | i will be shown In the table. 

* 

Along each C characteristic line (I « 1, 2,...) the data of bow shock, first 
point of cV after the bow shod?, the point of the formation of shock, and 
thcembedded shock whenever it is present will be tabulated. In case the 
same I of bow shock appears twice, that means the control In the step size 
is activated in the bow shock calculations (see remarks in regard to the step 
size control). The output format is as follows; 


Bow shock 


i * 


?' g" ir ?' s' h u- 

JMAX u v W V % | 

« sh 


Cpcm J - 
\ 


x" r" 0" H-r' 


sh 

'V 


% 

u v w y w s 


Siiock FORMED at J “ 
x" 


r" 0' IV V 

a a a ®a 


'V 

u 


a 


r\, 

v w 
a a 


KS = 


sir 


p(ks) 


C03 


2 


sh 


- i v- 


F SHOCK 


FORMED 


CASE{J 


KS » 


x' 

a 


0' 

n 

ti. 


°b 

u. 


% 


V 


H b 

a> 

v. 


w 


a i 


cob 


s r 

sh 


w, 


r 

$ . 

fih a 

r XMCAL 
L Xl‘ 

% 


“a 


C-* {' 

b b l XMDEFC 


y«,b 


sh, 


REMARK: Subscripts a, b for points ahead and behind the shock, 

XMDEFC and XMCAL are the difference of M i s computed from 
the shock equations and from the characteristic equations. 


EMSHOCK 


u 


°fc 

u, 


» 

H. 

c 

v 


v. 


r »d 


w 


“b 


w. 


KS » 

5 


sh. 


«a 


sh. 


»b 


r 

l 


sh. 


* 


sh, 


When the calculation reaches to the ground, all the data will be inter- 
polated. The output format is 


REACH ground. 


EMSHOCK 
, dpern 
cpSH 


x 0 M 


!\ 


fo <v 

U V w 


sh 


Ap 


5 sh 


(tfl) 


REMARK: Ap ~ p - in Tbs./sq. ft. 

Aq « pb for embedded shock only 


v- 


